{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "f9c7190b-f8e8-4df9-9f0b-36b1bed70542",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Loading required package: survival\n",
      "\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "Competing Risks Regression\n",
       "\n",
       "Call:\n",
       "crr(ftime = d$ftime, fstatus = d$fstatus, cov1 = X, failcode = 1, \n",
       "    cencode = 0)\n",
       "\n",
       "                          coef exp(coef) se(coef)     z p-value\n",
       "statin_baseline       -0.31863     0.727  0.18328 -1.74 8.2e-02\n",
       "age                    0.02930     1.030  0.00414  7.08 1.5e-12\n",
       "male                   0.22463     1.252  0.09269  2.42 1.5e-02\n",
       "bmi                    0.00607     1.006  0.00189  3.21 1.3e-03\n",
       "hb_glic                0.16402     1.178  0.02900  5.66 1.5e-08\n",
       "creatinine_first       0.34864     1.417  0.04914  7.09 1.3e-12\n",
       "hypertension_baseline  0.28320     1.327  0.10324  2.74 6.1e-03\n",
       "cvd_baseline           0.59807     1.819  0.10023  5.97 2.4e-09\n",
       "ckd_baseline           1.70113     5.480  0.14192 11.99 0.0e+00\n",
       "ace_arb_base           0.15811     1.171  0.15966  0.99 3.2e-01\n",
       "\n",
       "                      exp(coef) exp(-coef)  2.5% 97.5%\n",
       "statin_baseline           0.727      1.375 0.508  1.04\n",
       "age                       1.030      0.971 1.021  1.04\n",
       "male                      1.252      0.799 1.044  1.50\n",
       "bmi                       1.006      0.994 1.002  1.01\n",
       "hb_glic                   1.178      0.849 1.113  1.25\n",
       "creatinine_first          1.417      0.706 1.287  1.56\n",
       "hypertension_baseline     1.327      0.753 1.084  1.63\n",
       "cvd_baseline              1.819      0.550 1.494  2.21\n",
       "ckd_baseline              5.480      0.182 4.149  7.24\n",
       "ace_arb_base              1.171      0.854 0.857  1.60\n",
       "\n",
       "Num. cases = 36246\n",
       "Pseudo Log-likelihood = -5774 \n",
       "Pseudo likelihood ratio test = 1383  on 10 df,"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Tests:\n",
      "      stat          pv df\n",
      "1 7.517797 0.006109234  1\n",
      "2 4.594160 0.032081057  1\n",
      "Estimates and Variances:\n",
      "$est\n",
      "                               1000        2000        3000        4000\n",
      "No statin at baseline 1 0.017240979 0.028520282 0.042874830 0.056103138\n",
      "Statin at baseline 1    0.027733032 0.042945951          NA          NA\n",
      "No statin at baseline 2 0.003749436 0.003749436 0.003749436 0.003749436\n",
      "Statin at baseline 2    0.007646365 0.007646365          NA          NA\n",
      "\n",
      "$var\n",
      "                                1000         2000         3000         4000\n",
      "No statin at baseline 1 7.083441e-07 1.747125e-06 6.810305e-06 1.007876e-04\n",
      "Statin at baseline 1    2.171749e-05 7.231943e-05           NA           NA\n",
      "No statin at baseline 2 1.470002e-07 1.470002e-07 1.470002e-07 1.470002e-07\n",
      "Statin at baseline 2    5.676679e-06 5.676679e-06           NA           NA\n",
      "\n"
     ]
    }
   ],
   "source": [
    "library(cmprsk)\n",
    "\n",
    "d <- read.csv(\"finegray_input.csv\")\n",
    "\n",
    "# Matrix of covariable (remove intercept in matrix of model.matrix)\n",
    "X <- model.matrix(~ statin_baseline + age + male + bmi + hb_glic + creatinine_first +\n",
    "                    hypertension_baseline + cvd_baseline + ckd_baseline + ace_arb_base,\n",
    "                  data = d)[, -1]\n",
    "\n",
    "# Fine–Gray (evento 1 = diálise; 0 = censura; 2 = morte)\n",
    "fit_fg <- crr(ftime = d$ftime, fstatus = d$fstatus, cov1 = X, failcode = 1, cencode = 0)\n",
    "summary(fit_fg)   # sHR, IC95%, p\n",
    "\n",
    "# Gray's test + CIF por grupos de exposição no baseline\n",
    "g <- factor(d$statin_baseline, labels = c(\"No statin at baseline\", \"Statin at baseline\"))\n",
    "cif <- cuminc(ftime = d$ftime, fstatus = d$fstatus, group = g, cencode = 0)\n",
    "print(cif)        # inclui Gray's test\n",
    "# plot(cif, lty=c(1,2), col=1, xlab=\"Days since T2DM diagnosis\", ylab=\"Cumulative incidence of dialysis\")\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "f7f6a924-c03c-411d-a335-13336adc8af8",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "=== Fine–Gray: finegray_newuser_lag30_censor.csv ===\n",
      "Competing Risks Regression\n",
      "\n",
      "Call:\n",
      "crr(ftime = d$ftime, fstatus = d$fstatus, cov1 = X, failcode = 1, \n",
      "    cencode = 0)\n",
      "\n",
      "                         coef exp(coef) se(coef)     z p-value\n",
      "statin_newuser        -0.2240     0.799  0.12341 -1.82 6.9e-02\n",
      "age                    0.0336     1.034  0.00514  6.53 6.4e-11\n",
      "male                   0.2150     1.240  0.11469  1.87 6.1e-02\n",
      "bmi                    0.0071     1.007  0.00233  3.04 2.3e-03\n",
      "hb_glic                0.1568     1.170  0.03746  4.19 2.8e-05\n",
      "creatinine_first       0.3154     1.371  0.04579  6.89 5.6e-12\n",
      "hypertension_baseline  0.2279     1.256  0.12749  1.79 7.4e-02\n",
      "cvd_baseline           0.5089     1.663  0.13396  3.80 1.5e-04\n",
      "ckd_baseline           1.6882     5.410  0.19126  8.83 0.0e+00\n",
      "ace_arb_base           0.3063     1.358  0.22485  1.36 1.7e-01\n",
      "\n",
      "                      exp(coef) exp(-coef)  2.5% 97.5%\n",
      "statin_newuser            0.799      1.251 0.628  1.02\n",
      "age                       1.034      0.967 1.024  1.04\n",
      "male                      1.240      0.807 0.990  1.55\n",
      "bmi                       1.007      0.993 1.003  1.01\n",
      "hb_glic                   1.170      0.855 1.087  1.26\n",
      "creatinine_first          1.371      0.729 1.253  1.50\n",
      "hypertension_baseline     1.256      0.796 0.978  1.61\n",
      "cvd_baseline              1.663      0.601 1.279  2.16\n",
      "ckd_baseline              5.410      0.185 3.718  7.87\n",
      "ace_arb_base              1.358      0.736 0.874  2.11\n",
      "\n",
      "Num. cases = 34242\n",
      "Pseudo Log-likelihood = -3602 \n",
      "Pseudo likelihood ratio test = 665  on 10 df,\n",
      "\n",
      "=== Fine–Gray: finegray_newuser_lag90_censor.csv ===\n",
      "Competing Risks Regression\n",
      "\n",
      "Call:\n",
      "crr(ftime = d$ftime, fstatus = d$fstatus, cov1 = X, failcode = 1, \n",
      "    cencode = 0)\n",
      "\n",
      "                         coef exp(coef) se(coef)     z p-value\n",
      "statin_newuser        -0.3402     0.712  0.12626 -2.69 7.1e-03\n",
      "age                    0.0356     1.036  0.00525  6.78 1.2e-11\n",
      "male                   0.1987     1.220  0.11613  1.71 8.7e-02\n",
      "bmi                    0.0074     1.007  0.00230  3.21 1.3e-03\n",
      "hb_glic                0.1498     1.162  0.03895  3.85 1.2e-04\n",
      "creatinine_first       0.3171     1.373  0.04628  6.85 7.3e-12\n",
      "hypertension_baseline  0.2698     1.310  0.13020  2.07 3.8e-02\n",
      "cvd_baseline           0.4742     1.607  0.13779  3.44 5.8e-04\n",
      "ckd_baseline           1.6724     5.325  0.19622  8.52 0.0e+00\n",
      "ace_arb_base           0.3369     1.401  0.22549  1.49 1.4e-01\n",
      "\n",
      "                      exp(coef) exp(-coef)  2.5% 97.5%\n",
      "statin_newuser            0.712      1.405 0.556 0.911\n",
      "age                       1.036      0.965 1.026 1.047\n",
      "male                      1.220      0.820 0.972 1.532\n",
      "bmi                       1.007      0.993 1.003 1.012\n",
      "hb_glic                   1.162      0.861 1.076 1.254\n",
      "creatinine_first          1.373      0.728 1.254 1.504\n",
      "hypertension_baseline     1.310      0.764 1.015 1.690\n",
      "cvd_baseline              1.607      0.622 1.226 2.105\n",
      "ckd_baseline              5.325      0.188 3.625 7.823\n",
      "ace_arb_base              1.401      0.714 0.900 2.179\n",
      "\n",
      "Num. cases = 34242\n",
      "Pseudo Log-likelihood = -3464 \n",
      "Pseudo likelihood ratio test = 644  on 10 df,\n",
      "Tests:\n",
      "      stat           pv df\n",
      "1 18.69291 1.535528e-05  1\n",
      "2 10.18785 1.413687e-03  1\n",
      "Estimates and Variances:\n",
      "$est\n",
      "                             1000        2000        3000        4000\n",
      "Never users 1         0.010146945 0.014602588 0.022280362 0.022280362\n",
      "New users (lag 30d) 1 0.011467243 0.030409368 0.047623389 0.071850520\n",
      "Never users 2         0.004130607 0.004130607 0.004130607 0.004130607\n",
      "New users (lag 30d) 2 0.001727673 0.001727673 0.001727673 0.001727673\n",
      "\n",
      "$var\n",
      "                              1000         2000         3000         4000\n",
      "Never users 1         5.657544e-07 1.222755e-06 5.560925e-06 5.560925e-06\n",
      "New users (lag 30d) 1 1.872928e-06 6.953871e-06 2.165811e-05 3.153154e-04\n",
      "Never users 2         2.260075e-07 2.260075e-07 2.260075e-07 2.260075e-07\n",
      "New users (lag 30d) 2 2.499860e-07 2.499860e-07 2.499860e-07 2.499860e-07\n",
      "\n"
     ]
    }
   ],
   "source": [
    "# install.packages(\"cmprsk\", repos = \"https://cloud.r-project.org\")\n",
    "library(cmprsk)\n",
    "\n",
    "run_fg <- function(csv_path) {\n",
    "  d <- read.csv(csv_path)\n",
    "\n",
    "  # matriz de covariáveis (sem intercepto duplicado)\n",
    "  X <- model.matrix(~ statin_newuser + age + male + bmi + hb_glic + creatinine_first +\n",
    "                    hypertension_baseline + cvd_baseline + ckd_baseline + ace_arb_base,\n",
    "                    data = d)[, -1]\n",
    "\n",
    "  fit <- crr(ftime = d$ftime, fstatus = d$fstatus, cov1 = X, failcode = 1, cencode = 0)\n",
    "  cat(\"\\n=== Fine–Gray:\", csv_path, \"===\\n\")\n",
    "  print(summary(fit))  # sHR, IC95%, p\n",
    "  invisible(fit)\n",
    "}\n",
    "\n",
    "# rode para 30 e 90 dias\n",
    "fit30 <- run_fg(\"finegray_newuser_lag30_censor.csv\")\n",
    "fit90 <- run_fg(\"finegray_newuser_lag90_censor.csv\")\n",
    "\n",
    "# (opcional) CIF por grupos (new-user vs never) no baseline desta análise\n",
    "d30 <- read.csv(\"finegray_newuser_lag30_censor.csv\")\n",
    "g <- factor(d30$statin_newuser, labels = c(\"Never users\",\"New users (lag 30d)\"))\n",
    "cif <- cuminc(ftime = d30$ftime, fstatus = d30$fstatus, group = g, cencode = 0)\n",
    "print(cif)  # inclui Gray’s test\n",
    "# plot(cif, lty = c(1,2), col = 1, xlab = \"Days since T2DM diagnosis\", ylab = \"Cumulative incidence of dialysis\")\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a45c08d1-6c25-4562-a143-4bc41512511b",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "R",
   "language": "R",
   "name": "ir"
  },
  "language_info": {
   "codemirror_mode": "r",
   "file_extension": ".r",
   "mimetype": "text/x-r-source",
   "name": "R",
   "pygments_lexer": "r",
   "version": "4.5.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
